################################################################################
################################################################################
### Script file to reproduce the empirical results from the Strategic Ambiguity project (JEPOP version)
### Victor Shin and Laron K. Williams
###
### Created: 6-4-24
### Modified: 
###
### Dependencies: SW.RData
### Output: all figures in the manuscript and appendix
###
################################################################################


################################################################################
################################################################################

library(lme4)
library(ggplot2)
library(texreg)
library(marginaleffects)

################################################################################

# Clear the space
rm(list = ls())

################################################################################
################################## Data Management #############################
################################################################################

# Set up the working directory
#setwd("")

# Read in the RData which contains the two datasets
load("SW.RData")


################################################################################
################################## Summary Statistics ##########################
################################################################################

# Individual-level data
summary(data.ind[, c('votefor', 'dist', 'pdf', 'scoreA', 'brand', 'prtycore', 'like')])
sd(data.ind$votefor, na.rm = T)
sd(data.ind$dist, na.rm = T)
sd(data.ind$pdf, na.rm = T)
sd(data.ind$scoreA, na.rm = T)
sd(data.ind$brand, na.rm = T)
sd(data.ind$prtycore, na.rm = T)
sd(data.ind$like, na.rm = T)


# Party-level data
summary(data.party[, c('vote_cses', 'scoreA.x', 'brand', 'pdf', 'xtrm', 'enp', 'dist2')])
sd(data.party$vote_cses, na.rm = T)
sd(data.party$scoreA.x, na.rm = T)
sd(data.party$brand, na.rm = T)
sd(data.party$pdf, na.rm = T)
sd(data.party$xtrm, na.rm = T)
sd(data.party$enp, na.rm = T)
sd(data.party$dist2, na.rm = T)



################################################################################
################################## Models ######################################
################################################################################

### Manuscript
# Model 1
m1 <- lmer(dist ~ brand*pdf + prtycore + like + adifp + (1| ccses), data = data.ind)
summary(m1)  

# Model 2
m2 <- glmer(votefor ~ brand*pdf + prtycore + like + dist + (1|ccses), data = data.ind, family = binomial)
summary(m2)

texreg(list(m1, m2), stars = c(0.01, 0.05, 0.1), caption.above = T)

### Appendix
# Model 3
m3 <- glmer(votefor ~ brand*pdf + prtycore + like + adifp + (1|ccses), data = data.ind, family = binomial)
summary(m3)

texreg(list(m1, m2, m3), stars = c(0.01, 0.05, 0.1), caption.above = T)

# Model 4
m4 <- lmer(dist ~ scoreA*pdf + prtycore + like + adifp + (1| ccses), data = data.ind)
summary(m4)  

# Model 5
m5 <- glmer(votefor ~ scoreA*pdf + prtycore + like + dist + (1|ccses), data = data.ind, family = binomial)
summary(m5)

# Model 6
m6 <- glmer(votefor ~ scoreA*pdf + prtycore + like + adifp + (1|ccses), data = data.ind, family = binomial)
summary(m6)

texreg(list(m4, m5, m6), stars = c(0.01, 0.05, 0.1), caption.above = T)

# Model 7
m7 <- lmer(vote_cses ~ brand*pdf + dist2 + xtrm + enp + lagvote + (1|dyad), data = data.party)
summary(m7)

# Model 8
m8 <- lmer(vote_cses ~ scoreA.x*pdf + dist2 + xtrm + enp + lagvote + (1|dyad), data = data.party)
summary(m8)

texreg(list(m7, m8), stars = c(0.01, 0.05, 0.1), caption.above = T)







################################################################################
################################## Figures #####################################
################################################################################

### Manuscript
# Model 1: Marginal effect of ambiguity (brand dispersion) on perceived distance across ideological overlap
pdf("Manuscript/Figures/m1_plot_a.pdf")
m1.plot.a <- plot_comparisons(m1,
                             variable = list(brand="sd"),
                             condition = "pdf", re.form = NA) +
  labs(x = "Ideological Overlap", y = "Perceived Distance") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=pdf), data=data.ind, sides="b") +
  theme_classic()
m1.plot.a
dev.off()


# Model 2: Marginal effect of ambiguity (brand dispersion) on vote choice across ideological overlap
pdf("Manuscript/Figures/m2_plot_a.pdf")
m2.plot.a <- plot_comparisons(m2,
                             variable = list(brand="sd"),
                             condition = "pdf", re.form = NA) +
  labs(x = "Ideological Overlap", y = "Pr(Party Vote)") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=pdf), data=data.ind, sides="b") +
  theme_classic()
m2.plot.a
dev.off()


# Model 2: Marginal effect of ideological overlap on vote choice across ambiguity (brand dispersion)
pdf("Manuscript/Figures/m2_plot_b.pdf")
m2.plot.b <- plot_comparisons(m2,
                              variable = list(pdf="sd"),
                              condition = "brand", re.form = NA) +
  labs(x = "Ambiguity", y = "Pr(Party Vote)") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=brand), data=data.ind, sides="b") +
  theme_classic()
m2.plot.b
dev.off()



### Appendix
# Model 3: Marginal effect of ambiguity (brand dispersion) on vote choice across ideological overlap
pdf("Appendix/Figures/m3_plot_a.pdf")
m3.plot.a <- plot_comparisons(m3,
                              variable = list(brand="sd"),
                              condition = "pdf", re.form = NA) +
  labs(x = "Ideological Overlap", y = "Pr(Party Vote)") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=pdf), data=data.ind, sides="b") +
  theme_classic()
m3.plot.a
dev.off()


# Model 3: Marginal effect of ideological overlap on vote choice across ambiguity (brand dispersion)
pdf("Appendix/Figures/m3_plot_b.pdf")
m3.plot.b <- plot_comparisons(m3,
                              variable = list(pdf="sd"),
                              condition = "brand", re.form = NA) +
  labs(x = "Ambiguity", y = "Pr(Party Vote)") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=brand), data=data.ind, sides="b") +
  theme_classic()
m3.plot.b
dev.off()



# Model 4: Marginal effect of ambiguity on perceived distance across ideological overlap
pdf("Appendix/Figures/m4_plot_a.pdf")
m4.plot.a <- plot_comparisons(m4,
                             variable = list(scoreA = "sd"),
                             condition = "pdf", re.form = NA) +
  labs(x = "Ideological Overlap", y = "Perceived Distance") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=pdf), data=data.ind, sides="b") +
  theme_classic() 
m4.plot.a
dev.off()



# Model 5: Marginal effect of ambiguity on vote choice across ideological overlap
pdf("Appendix/Figures/m5_plot_a.pdf")
m5.plot.a <- plot_comparisons(m5,
                             variable = list(scoreA="sd"),
                             condition = "pdf", re.form = NA) +
  labs(x = "Ideological Overlap", y = "Pr(Party Vote)") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=pdf), data=data.ind, sides="b") +
  theme_classic()
m5.plot.a
dev.off()


# Model 5: Marginal effect of ideological overlap on vote choice across ambiguity
pdf("Appendix/Figures/m5_plot_b.pdf")
m5.plot.b <- plot_comparisons(m5,
                             variable = list(pdf="sd"),
                             condition = "scoreA", re.form = NA) +
  labs(x = "Ambiguity", y = "Pr(Party Vote)") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=scoreA), data=data.ind, sides="b") +
  theme_classic()
m5.plot.b
dev.off()


# Model 6: Marginal effect of ambiguity on vote choice across ideological overlap
pdf("Appendix/Figures/m6_plot_a.pdf")
m6.plot.a <- plot_comparisons(m6,
                             variable = list(scoreA="sd"),
                             condition = "pdf", re.form = NA) +
  labs(x = "Ideological Overlap", y = "Pr(Party Vote)") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=pdf), data=data.ind, sides="b") +
  theme_classic()
m6.plot.a
dev.off()


# Model 6: Marginal effect of ideological overlap on vote choice across ambiguity
pdf("Appendix/Figures/m6_plot_b.pdf")
m6.plot.b <- plot_comparisons(m6,
                             variable = list(pdf="sd"),
                             condition = "scoreA", re.form = NA) +
  labs(x = "Ambiguity", y = "Pr(Party Vote)") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=scoreA), data=data.ind, sides="b") +
  theme_classic()
m6.plot.b
dev.off()



# Model 7: Marginal effect of ambiguity (brand dispersion) on vote share across ideological overlap
pdf("Appendix/Figures/m7_plot_a.pdf")
m7.plot.a <- plot_comparisons(m7, 
                              variables = list(brand = "sd"), 
                              condition = "pdf", re.form = NA) +
  labs(x = "Ideological Overlap", y = "Vote Share") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=pdf), data=data.party, sides="b") +
  xlim(0.1, 1) +
  theme_classic()
m7.plot.a
dev.off()


# Model 8: Marginal effect of ambiguity (disagreement) on vote share across ideological overlap
pdf("Appendix/Figures/m8_plot_a.pdf")
m8.plot.a <- plot_comparisons(m8, 
                              variables = list(scoreA.x = "sd"), 
                              condition = "pdf", re.form = NA) +
  labs(x = "Ideological Overlap", y = "Vote Share") +
  geom_hline(yintercept=0,linetype=2) +
  geom_rug(aes(x=pdf), data=data.party, sides="b") +
  xlim(0.1, 1) +
  theme_classic()
m8.plot.a
dev.off()



################################################################################
### Quantities of interest: set up the scenario values 
################################################################################

# Set up the scenarios: get the standard deviations and means
s0.scoreA <- mean(data7_adv$scoreA, na.rm = T)
sd.scoreA <- sd(data7_adv$scoreA, na.rm = T)
s1.scoreA <- s0.scoreA + sd.scoreA
min.scoreA <- min(data7_adv$scoreA, na.rm = T)
max.scoreA <- max(data7_adv$scoreA, na.rm = T)
s0.scoreA
sd.scoreA
s1.scoreA

s0.brand <- mean(data7_adv$scoreA, na.rm = T)
sd.brand <- sd(data7_adv$scoreA, na.rm = T)
s1.brand <- s0.brand + sd.brand
min.brand <- min(data7_adv$scoreA, na.rm = T)
max.brand <- max(data7_adv$scoreA, na.rm = T)
s0.brand
sd.brand
s1.brand

s0.pdf <- mean(data7_adv$pdf, na.rm = T)
sd.pdf <- sd(data7_adv$pdf, na.rm = T)
s1.pdf <- s0.pdf + sd.pdf
min.pdf <- min(data7_adv$pdf, na.rm = T)
max.pdf <- max(data7_adv$pdf, na.rm = T)
s0.pdf
sd.pdf
s1.pdf

s0.prtycore <- 0
s1.prtycore <- 1

s0.like <- mean(data7_adv$like, na.rm = T)
sd.like <- sd(data7_adv$like, na.rm = T)
s1.like <- s0.like + sd.like
s0.like
sd.like
s1.like

s0.adifp <- mean(data7_adv$adifp, na.rm = T)
sd.adifp <- sd(data7_adv$adifp, na.rm = T)
s1.adifp <- s0.adifp + sd.adifp
s0.adifp
sd.adifp
s1.adifp

s0.dist <- mean(data7_adv$dist, na.rm = T)
sd.dist <- sd(data7_adv$dist, na.rm = T)
s1.dist <- s0.dist + sd.dist
s0.dist
sd.dist
s1.dist

################################################################################
### Quantities of interest
################################################################################

# Model 1
m1 <- lmer(dist ~ brand*pdf + prtycore + like + adifp + (1| ccses), data = data7_adv)
summary(m1)

m1.sim <- arm::sim(m1, 1000)

# What is the effect of partisanship?
effect <- (m1.sim@fixef[ , "prtycore"]) * 1

m1.prtycore <- data.frame(
  me = quantile(effect, 0.5),
  lo.95 = quantile(effect, 0.025),
  hi.95 = quantile(effect, 0.975),
  lo.90 = quantile(effect, 0.05),
  hi.90 = quantile(effect, 0.95)
)
m1.prtycore

# What is the effect of likeability?
effect <- (m1.sim@fixef[ , "like"]) * sd.like

m1.like <- data.frame(
  me = quantile(effect, 0.5),
  lo.95 = quantile(effect, 0.025),
  hi.95 = quantile(effect, 0.975),
  lo.90 = quantile(effect, 0.05),
  hi.90 = quantile(effect, 0.95)
)
m1.like

# What is the effect of actual distance?
effect <- (m1.sim@fixef[ , "adifp"]) * sd.adifp

m1.adifp <- data.frame(
  me = quantile(effect, 0.5),
  lo.95 = quantile(effect, 0.025),
  hi.95 = quantile(effect, 0.975),
  lo.90 = quantile(effect, 0.05),
  hi.90 = quantile(effect, 0.95)
)
m1.adifp




# Model 2
m2 <- glmer(votefor ~ brand*pdf + prtycore + dist + like + (1|ccses), data = data7_adv, family = binomial)
summary(m2)
m2.sim <- arm::sim(m2, 1000)

# Average probability across the sample
avg.prob <- apply(fitted(m2.sim, m2), 1, function(x) quantile(x, 0.5))
mean(avg.prob)

# All the other variables at their means (or modes)

# What is the effect of partisanship?
s0 <- plogis(
        m2.sim@fixef[, "(Intercept)"] + m2.sim@fixef[, "brand"] * s0.brand + 
        m2.sim@fixef[, "pdf"] * s0.pdf + m2.sim@fixef[, "prtycore"] * 0 +
        m2.sim@fixef[, "dist"] * s0.dist + m2.sim@fixef[, "like"] * s0.like + 
        m2.sim@fixef[, "brand:pdf"] * s0.brand * s0.pdf
      )

s1 <- plogis(
  m2.sim@fixef[, "(Intercept)"] + m2.sim@fixef[, "brand"] * s0.brand + 
    m2.sim@fixef[, "pdf"] * s0.pdf + m2.sim@fixef[, "prtycore"] * 1 +
    m2.sim@fixef[, "dist"] * s0.dist + m2.sim@fixef[, "like"] * s0.like + 
    m2.sim@fixef[, "brand:pdf"] * s0.brand * s0.pdf
)

effect <- s1 - s0

m2.prtycore <- data.frame(
  me = quantile(effect, 0.5),
  lo.95 = quantile(effect, 0.025),
  hi.95 = quantile(effect, 0.975),
  lo.90 = quantile(effect, 0.05),
  hi.90 = quantile(effect, 0.95)
)
m2.prtycore


# What is the effect of likeability?
s0 <- plogis(
  m2.sim@fixef[, "(Intercept)"] + m2.sim@fixef[, "brand"] * s0.brand + 
    m2.sim@fixef[, "pdf"] * s0.pdf + m2.sim@fixef[, "prtycore"] * 0 +
    m2.sim@fixef[, "dist"] * s0.dist + m2.sim@fixef[, "like"] * s0.like + 
    m2.sim@fixef[, "brand:pdf"] * s0.brand * s0.pdf
)

s1 <- plogis(
  m2.sim@fixef[, "(Intercept)"] + m2.sim@fixef[, "brand"] * s0.brand + 
    m2.sim@fixef[, "pdf"] * s0.pdf + m2.sim@fixef[, "prtycore"] * 0 +
    m2.sim@fixef[, "dist"] * s0.dist + m2.sim@fixef[, "like"] * s1.like + 
    m2.sim@fixef[, "brand:pdf"] * s0.brand * s0.pdf
)

effect <- s1 - s0

m2.like <- data.frame(
  me = quantile(effect, 0.5),
  lo.95 = quantile(effect, 0.025),
  hi.95 = quantile(effect, 0.975),
  lo.90 = quantile(effect, 0.05),
  hi.90 = quantile(effect, 0.95)
)
m2.like

# What is the effect of perceived distance?
s0 <- plogis(
  m2.sim@fixef[, "(Intercept)"] + m2.sim@fixef[, "brand"] * s0.brand + 
    m2.sim@fixef[, "pdf"] * s0.pdf + m2.sim@fixef[, "prtycore"] * 0 +
    m2.sim@fixef[, "dist"] * s0.dist + m2.sim@fixef[, "like"] * s0.like + 
    m2.sim@fixef[, "brand:pdf"] * s0.brand * s0.pdf
)

s1 <- plogis(
  m2.sim@fixef[, "(Intercept)"] + m2.sim@fixef[, "brand"] * s0.brand + 
    m2.sim@fixef[, "pdf"] * s0.pdf + m2.sim@fixef[, "prtycore"] * 0 +
    m2.sim@fixef[, "dist"] * s1.dist + m2.sim@fixef[, "like"] * s0.like + 
    m2.sim@fixef[, "brand:pdf"] * s0.brand * s0.pdf
)

effect <- s1 - s0

m2.dist <- data.frame(
  me = quantile(effect, 0.5),
  lo.95 = quantile(effect, 0.025),
  hi.95 = quantile(effect, 0.975),
  lo.90 = quantile(effect, 0.05),
  hi.90 = quantile(effect, 0.95)
)
m2.dist


### At what values of ideological overlap does the effect of ambiguity flip from negative to positive (Model 1)?
m1.plot.a <- plot_comparisons(m1, variable = list(brand="sd"), condition = "pdf", re.form = NA, draw = FALSE) 
m1.plot.a

# flips at pdf = 0.595 (30th percentile)
quantile(data.ind$pdf, na.rm = T)
quantile(data.ind$pdf, c(.31), na.rm = T)

### At what values of ideological overlap does the effect of ambiguity flip from positive to negative (Model 2)?
m2.plot.a <- plot_comparisons(m2, variable = list(brand="sd"), condition = "pdf", re.form = NA, draw = FALSE)
m2.plot.a

# flips at pdf = 0.74 (61st percentile)
quantile(data.ind$pdf, na.rm = T)
quantile(data.ind$pdf, c(.61), na.rm = T)

################################################################################
### Justification for Multi-Level Analysis
################################################################################

lme0 <- lmer(vote_cses ~ (1 | dyad), data = data.party)
summary(lme0)

icc <- 29.49 / (29.49 + 81.15)
icc   # 26.7%

# ANOVA test
lme3a <- lmer(vote_cses ~ lagvote + scoreA.x*pdf + xtrm + dist2 + enp + (1 | dyad), data = data.party, REML=FALSE)
lme3b <- lmer(vote_cses ~ lagvote + scoreA.x*pdf + xtrm + dist2 + enp + (scoreA.x | dyad), data = data.party, REML=FALSE)
lme3c <- lmer(vote_cses ~ lagvote + scoreA.x*pdf + xtrm + dist2 + enp + (scoreA.x + xtrm | dyad), data = data.party, REML=FALSE)
lme3d <- lmer(vote_cses ~ lagvote + scoreA.x*pdf + xtrm + dist2 + enp + (scoreA.x + xtrm + dist2 | dyad), data = data.party, REML=FALSE)

mod3 <- anova(lme3a,lme3b,lme3c, lme3d)
mod3
